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Wave  Data  Processing  and  Analysis 
Tools  for  Developing  Wave  Boundary 
Forcing  for  CMS-Wave  and  GenCade 
Numerical  Models:  Part  1 

by  Kenneth  J.  Connell  and  Rusty  Permenter 


PURPOSE:  This  Coastal  and  Hydraulics  Engineering  Technical  Note  (CHETN)  presents  a 
collection  of  wave  data  processing  and  analysis  tools  for  assimilating  measured  wave  data  into 
wave  climatology  statistics  and  representative  data  subsets.  These  representative  data  subsets  can 
be  applied  to  force  the  offshore  boundary  of  numerical  wave  models  and  also  provide  general  wave 
climatology  analysis.  This  CHETN  focuses  on  applications  employing  the  numerical  models 
CMS-Wave  (Lin  et  al.  2008)  and  GenCade  (Frey  et  al.  2012);  however,  the  toolkit  has  been 
developed  with  substantial  flexibility  such  that  it  can  be  readily  modified  for  other  models  or 
analyses.  The  CHETN  first  describes  the  data  processing  tools  and  then  presents  an  example 
application.  The  general  work-flow  process,  specific  inputs,  execution  steps,  and  resulting  outputs 
are  described  through  the  example  application.  The  toolkit  and  the  example  dataset  discussed 
herein  can  be  downloaded  from  http://cirp.mace. army, mil/products/?tab=4. 

DESCRIPTION  OF  TOOLKIT:  The  wave  data  processing  toolkit  was  developed  to  operate  in 
the  MATLAB®  environment  (2007  or  later)  with  the  exception  of  the  final  step  of  CMS-Wave- 
to-GenCade  map  file  conversion  which  was  developed  in  Fortran.  The  toolkit  consists  of  three 
MATLAB  m-files  (* *.m)  and  one  Windows  Fortran  executable,  which  is  discussed  in  a  second 
technical  note  (Permenter  et  al.  2013): 

•  JP_BinWaves.m:  This  is  the  main  code  and  includes  input  options,  calls  to  the  other 
functions,  the  main  probability  histogram  binning  routine,  tabular  output  commands,  and 
graphical  data  output  plots. 

•  Linwavetheory.m:  This  function  calculates  the  linear  wave  theory  transformation  and 
output  and  calls  wavdisp  eq.m  function. 

•  Wavdispeq.m:  This  function  calculates  the  dispersion  equation  using  the  Newton-Raphson 
method  for  convergence. 

The  toolkit  process  assumes  the  user  has  downloaded  a  long-term,  time  series  wave  dataset  in  a 
text  file  (*.txt  or  *.dat)  in  meters  and  meteorological  wave  direction  convention.  These  data  were 
then  imported  into  MATLAB  and  processed  and  analyzed  during  the  execution  of  JP  WaveBin.m. 
The  tool  accepts  input  through  a  series  of  command  prompt  questions  in  order  to  determine  the 
desired  wave  binning  strategy. 

The  code  filters  the  data  and  then  bins  the  data  to  develop  a  3 -way  (significant  wave  height  \H/\  vs. 
wave  period  [7]  vs.  wave  directions  \Dir])  joint  probability  ( JP ).  The  code  also  allows  for  seasonal 


Approved  for  public  release;  distribution  is  unlimited. 


US  Army  Corps 
of  Engineers® 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

SEP  2013  2' REPORT  TYPE 

3.  DATES  COVERED 

00-00-2013  to  00-00-2013 

4.  TITLE  AND  SUBTITLE 

Wave  Data  Processing  and  Analysis  Tools  for  Developing  Wave 

Boundary  Forcing  for  CMS- Wave  and  GenCade  Numerical  Models:  Part 

1 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS (ES) 

US  Army  Engineer  Research  and  Development  Center  (ERDC), Coastal 
and  Hydraulics  Laboratory  (CHL), Vicksburg, MS, 39180 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS  (ES) 

10.  SPONSOR/MONITOR’S  ACRONYM(S) 

11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF:  17.  LIMITATION  OF 

ARSTRATT 

1 8 .  NUMBER  1 9a.  NAME  OF 

OF  PAGES  RESPONSIBLE  PERSON 

a.  REPORT  b.  ABSTRACT  c.  THIS  PAGE  Same  3S 

unclassified  unclassified  unclassified  Report  (SAR) 

10 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


ERDC/CHL  CHETN-IV-97 
September  2013 

binning  of  the  time  series  by  evenly  dividing  each  year  into  a  given  number  of  time  bins.  Tables 
for  two-way  JP  are  written  as  output  for  comparison  of  Hs  vs.  T,  Hs  vs.  Dir,  and  T  vs.  Dir.  The 
two-way  JP  tables  may  be  used  for  developing  wave  roses  using  the  bin  definitions,  the  bin  data 
(for  Hs,  T,  and  Dir),  and  the  percent  occurrence  (or  the  number  of  occurrences). 

The  user  defines  the  input  bins  for  Hs,  T,  and  Dir.  Similarly,  the  user  may  define  maximum  and 
minimum  values  for  wave  parameters  (e.g.,  only  include  events  with  Hs  above  3  m  or  disregard  all 
wave  events  with  T  less  than  4  sec;  or  Hs  less  than  0.25  m,  etc.)  to  remove  outliers  or  flags  in  the 
input  dataset.  The  code  generates  an  output  file  containing  the  percent  occurrence  for  each  time 
division  for  each  binned  wave  condition.  The  code  divides  each  year  of  the  time  series  into  equally 
spaced  intervals  to  allow  for  characterization  of  seasonal  variability  in  the  input  wave  time  series. 
The  user  specifies  the  number  of  seasonal  bins  used  (1  for  no  binning,  12  for  monthly).  The 
number  of  occurrences  for  each  height,  period,  and  direction  combination  is  calculated  for  each 
seasonal  bin  for  the  entire  length  of  the  input  dataset  (an  event  in  February  1991  will  have  the  same 
effects  as  an  equivalent  event  in  February  1995).  Each  row  in  the  output  file  represents  one  binned 
condition,  while  the  columns  represent  the  seasonal  bins.  In  addition  to  the  JP  table  output,  output 
also  includes  general  wave  statistics  (Mean  significant  wave  height,  Mean  peak  Period,  Standard 
Deviation,  Maximum  significant  wave  height,  Minimum  significant  wave  height,  and  Period  and 
Direction  associated  with  maximum  significant  wave  height)  and  a  file  containing  the  wave  index, 
height,  period,  direction,  number  of  occurrences,  and  percent  occurrences  for  each  binned  wave 
condition  without  seasonal  binning.  The  index  for  each  wave  condition  is  represented  by  the 
following  formula: 

Fndmx  nf  Event.  =  Index  nf  H  Rin  *  1  fiOOfl  +  Index,  nf  T  Bin.  *  1 00  +  In.d.ex  nf  D  Bin 

The  user  must  use  the  SMS  interface  to  create  a  *.eng  file  in  order  to  run  CMS-Wave.  The  table 
listed  in  the  WaveBin_Table_Over_MinCount_Threshold.dat  output  file  must  be  copied  into  the 
spectra  generator  within  CMS-Wave.  This  process  will  be  detailed  in  CHETN-IV-98.  The  *.eng 
file  supplies  the  offshore  wave  spectra  for  the  CMS-Wave  runs  (Lin  et  al.  2008).  The  CMS2MAP 
executable  takes  the  CMS-Wave  output  files  and  produces  a  *.map  or  *.gen  file  for  input  into 
GENCADE.  Features  of  the  wave  data  processing  toolkit  are  presented  below  in  an  example 
application  following  step-by-step  instructions  with  comments  to  guide  users. 

EXAMPLE:  Wave  Information  for  Ship  Island,  Mississippi.  This  example  demonstrates 
use  of  the  wave  data  processing  toolkit  to  transform  offshore  wave  data  to  the  model  boundary, 
calculate  the  wave  data  statistics,  and  produce  the  input  data  files  for  wave  models.  For  this 
example  the  user  has  saved  a  20-year  wave  time  series  from  WIS  station  144  at  a  depth  of  16  m. 
The  user  executes  JP_BinWaves.m  in  order  to  transform  the  waves  to  the  offshore  boundary  of  a 
CMS-Wave  grid  and  bin  the  data.  The  user  then  wishes  to  develop  boundary  wave  forcing  which 
characterizes  the  20-year  time  series  into  a  manageable  set  of  steady-state  forcing  parameters, 
which  will  be  converted  to  input  wave  energy  spectra.  Note  that  the  code  is  developed  for  metric 
units  with  wave  angles  in  degrees  and  conforming  to  meteorological  direction  convention.  The 
toolkit  files  and  Ship  Island  example  datasets  may  be  downloaded  from  http://cirp. usace. army. mil/ 
products/?  tab=4. 
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Input  of  Wave  Information. 

Step  1.  Save  the  toolkit  MATLAB  m- files  and  offshore  wave  data  into  the  same  working 
directory. 

Step  2.  Launch  MATLAB  and  execute  JP_BinWaves.m  by  typing  “JP  BinWaves”  at  the 
command  prompt.  The  following  is  the  command  prompt  when  the  user  executes  the  main  script: 

>  >  >JP_Bin  Waves 

Step  3.  Executing  the  MATLAB  script  initially  opens  a  dialogue  window  that  gives  the  user  the 
opportunity  to  browse  for  and  select  the  wave  input  file.  The  file  should  be  in  either  an  ASCII  *.txt 
or  *.dat  format.  When  it  is  selected,  the  user  selects  <OK>  and  the  data  file  is  loaded  into 
MATLAB. 

Step  4.  The  user  is  prompted  to  input  the  file  column  delimiter.  In  this  example,  the  file  was 
space  delimited: 

Enter  wave  data  file  column  delimiter  [space]: 

<  >  for  space,  <,  >  for  comma,  < \t>  for  tab,  <enter_delimiter_character(s)_here>  for 
custom: 


Step  5.  The  user  is  prompted  to  input  the  number  of  header  rows  (if  any).  In  this  example,  there 
was  one  header  row: 

Enter  number  of  header  rows  in  wave  data  file  [default  =  ()]: 

Setting  wavedata_delimiter  to  default  space  delimited. 

Number  of  Header  Rows  in  Wave  Data  File: 

1 

Step  6.  The  user  is  prompted  to  input  the  column  numbers  associated  with  the  height,  period,  and 
direction  as  well  as  the  data  input  type  and  number  of  seasonal  bins.  The  seasonal  binning  is 
performed  by  evenly  dividing  the  year  into  the  desired  number  of  bins  and  determining  the  percent 
occurrence  of  each  binned  event  in  each  seasonal  bin.  The  data  input  type  is  used  to  determine 
which  columns  in  the  data  file  signify  the  date  stamps  (YYYYMMDDHHMMSS  is  the  first 
column  in  WIS).  This  step  provides  flexibility  to  change  format  of  the  input  data.  In  this  example, 
columns  10,  11,  and  16  represent  Hs,  T,  and  Dir,  respectively,  a  WIS  data  file  is  utilized,  and  12 
seasonal  bins  are  specified: 

Enter  column  number  for  Significant  Wave  Height  (H,): 

»>  10 

Enter  column  number  for  Wave  Period  (T): 

»>  11 

Enter  column  number  for  Wave  Direction  (Dir): 
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»>  16 

Enter  input  type  (1  for  NDBC  buoy,  2  for  WIS  Data  file): 

Enter  number  of  seasonal  periods  for  binning: 

»>  12 

Step  7.  At  this  point,  the  program  reads  the  data  file  into  MATLAB. 

Transformation  to  a  CMS-Wave  Offshore  Boundary. 

Step  1.  First,  the  user  is  prompted  to  make  a  decision  on  whether  or  not  to  apply  a  linear  wave 
transformation  from  a  distant  offshore  data  point  to  the  shallower  and  further  inshore  wave  model 
boundary.  In  cases  where  the  user  wishes  to  simply  process  and  analyze  the  data  directly  from  the 
source,  wave  transformation  may  not  be  necessary  or  desired.  However,  if  the  user  wishes  to 
translate  waves  from  a  buoy  (or  hindcast  station)  that  is  not  located  at  the  offshore  boundary  of  the 
wave  model  grid,  it  is  recommended  that  a  linear  transformation  be  applied  prior  to  continuing 
with  the  analysis.  In  this  example,  the  user  requests  linear  wave  transformation  to  shallower  water: 

Perform  Linear  Wave  Transformation  Routine  ( Snell's  Law)  <Y/[N]>? 

Step  2.  The  user  is  prompted  to  input  the  depth  of  the  offshore  input  wave  information  and  the 
desired  target  depth  to  transform  the  waves.  In  this  example,  the  offshore  depth  at  WIS  144  is  16 
m  and  the  target  depth  of  the  offshore  boundary  of  the  CMS-Wave  grid  is  12  m: 

Enter  input  wave  information  depth  (m): 

»>  16 

Enter  target  wave  station  depth  to  transform  waves  to  (m): 

»>12 

Step  3.  The  user  is  prompted  to  convert  the  wave  angle  from  meteorological  to  shore  normal.  In 
order  to  use  the  linear  wave  transformation  code,  the  waves  must  be  between  -90  and  90  degrees 
relative  to  shore  normal  (Figure  1). 

Step  4.  The  user  is  prompted  to  enter  the  wave  angle  of  the  shore-normal  reference: 

Enter  shore  normal  reference  angle  (e.g.,  270  =  west  facing  coast;  0  =  north  facing 
coast;  90  =  east  facing  coast;  180  =  south  facing  coast): 

»>  161.3 

Step  5.  The  user  is  prompted  to  enter  cutoff  values  for  both  the  wave  height  and  wave  period  to 
discard  outliers  and  flags  for  missing  data.  Any  direction  values  greater  than  90  or  less  than  -90 
will  be  automatically  discarded  at  this  point. 

Enter  threshold  of  maximum  wave  height: 

(all  H  values  greater  than  this  threshold  will  be  reset  to  NaN) 

»>  25 
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Figure  1 .  Shore  Normal  Wave  Convention 


Wave  Heights  >  25  have  been  reset  to  NaN! 

Enter  threshold  of  maximum  wave  period: 

(all  T  values  greater  than  this  threshold  will  be  reset  to  NaN) 

»>  20 

Wave  Period  >20  have  been  reset  to  NaN! 

Wave  Angles  >  90  have  been  reset  to  NaN! 

Wave  Angles  <-90  have  been  reset  to  NaN! 

Step  6.  After  the  angles  are  converted  to  shore-normal  reference  system,  the  linear  wave 
transformation  routine  is  executed.  This  is  where  the  call  to  linwavetheory.m  is  executed: 

Continue  with  calculations?  <Y/[N]> 

>>>  y 

Step  7.  Next,  the  call  to  wavdispeq.m  is  automatically  executed  from  within  the  linwavetheory.m 
function  in  order  to  calculate  the  wave  number  ( k )  at  the  target  station  depth. 

Step  8.  Finally,  the  wave  transformation  is  complete  and  the  routine  automatically  proceeds  to  the 
next  step:  histogram  bin  sorting  of  the  waves  and  calculating  mean  conditions  of  the  population  in 
each  bin. 

Time  Series  Conversion  to  Representative  Wave  Forcing. 

Step  1.  The  user  is  now  prompted  to  select  which  method  to  apply  to  define  the  histogram  binning 
scheme.  The  first  option  is  to  input  band  widths  and  auto-bin  between  minimum  and  maximum 
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data  values.  If  the  range  of  data  does  not  come  to  an  even  number  of  bins,  the  last  bin  contains  the 
data  from  the  end  of  the  previous  bin  to  the  maximum  data  value.  For  the  second  method  the  user 
inputs  a  bracketed  series  of  the  bin  edges  for  H,  T,  and  Direction  (e.g.,  [-90  -60  -30  -15  15  30  60 
90]).  In  this  example,  the  default  method  (method  1)  is  selected  such  that  the  user  may  define  the 
bin  width  for  auto-binning  between  min  and  max  data  values: 

Enter  input  method  to  define  bin  edges  (Input  <[1],  2>):1 

Input  band  widths  for  auto-binning  between  min  and  max  data  values;  Enter  1: 

Input  explicit  array  of  bin  edges;  Enter  2: 

In  this  example,  significant  wave  height  is  binned  in  0.5-m  increments;  wave  period  bins  in  2  sec 
increments;  and  wave-angle  bins  ranging  from  -  10-deg  increments. 

Enter  user-specified  wave  height  band  width  (e.g.,  0.5)  :  .5 
Enter  user-specified  wave  period  band  width  (e.g.,  2)  :  2 
Enter  user-specified  wave  angle  band  width  (e.g.,  10)  :  10 

Step  2.  The  histogram  bin  scheme  is  echoed  back  to  the  screen  so  that  the  user  may  review  the 
calculated  bin  scheme  based  on  the  user-defined  bin  bounds. 

Step  3.  Finally,  the  main  loop  is  executed  and  the  data  were  binned  into  each  user-defined 
histogram  bin  by  frequency  of  occurrence  of  each  permutation. 

Step  4.  An  option  allows  the  user  to  reduce  the  number  of  permutations  for  modeling  by 
ignoring  permutations  that  fall  below  a  minimum  count  threshold.  In  this  example,  the  user 
decides  to  include  every  permutation  occurring: 

Enter  minimum  count  (occurrence)  threshold  to  consider  for  output  table  <[1]>: 

Step  5.  After  the  main  binning  loop  is  complete,  a  series  of  ASCII  output  tables  reporting  results 
of  the  data  processing  executed  by  the  script  were  written  to  output  files  in  the  working  directory 
and  include  the  following  seven  files: 

•  WaveBin_table.dat:  Wave  bin  3-way  joint  probability  permutation  table 

•  WaveBin_  table_Over_MinCount_Threshold.dat:  Wave  bin  permutation  table  eliminating 
values  less  than  the  minimum  count  threshold 

•  WaveStats_table.dat:  Wave  data  bulk  statistics  table 

•  JP_HvsT_table.dat:  Significant  wave  height  vs.  period  2-way  joint  probability  table 

•  JP_DvsH_table.dat:  Wave  direction  vs.  significant  wave  height  2-way  joint  probability  table 

•  JP_TvsD_table.dat:  Wave  period  vs.  direction  2-way  joint  probability  table 

•  JP_MultipleTimebin_Probability.dat:  Wave  Permutation  table  with  Seasonal  Bins 

Step  6.  Finally,  a  series  of  figures  representing  the  data  processing  executed  by  the  script  are 
automatically  plotted  as  output.  These  figures  include  stacked  histogram  plots  that  show  the 
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probability  distribution  of  each  two-way  joint  probability  result  for  the  binned  significant  wave 
height,  period,  and  direction;  3-dimensional  (3D)  surface  histogram  plots  that  show  the  probability 
distribution  of  each  2-way  result  for  Hs  vs.  T  (demonstrated  for  the  example  case  in  Figure  2), 
wave  angle  to  Hs,  and  T  to  wave  angle;  and,  a  2D  probability  plot  of  Hs  vs.  T. 


Figure  2.  Percent  Occurrence  for  Wave  Height  vs.  Wave  Period 


At  this  stage,  the  user  is  prompted  on  whether  to  plot  the  polar  axis  of  the  wave  rose  as 
referencing  frequency  in  counts  or  frequency  by  percent  (%)  occurrence.  Percent  occurrence  is 
the  default  option: 

Enter  1  to  plot  frequency  in  counts  or  Enter  2  to  plot  %  Occurrence,  [2] : 

Plotting  frequency  in  %  Occurrence 
You  entered:  2 

Next,  a  wave  rose  is  developed  and  plotted  as  shown  in  Figure  3.  It  is  important  to  note  that  the 
wave  directions  are  referenced  to  shore  normal  and  therefore  fall  between  -90  deg  and  +90  deg. 
The  polar  axis  may  be  rotated  to  any  angle  with  minimal  effort  within  MATLAB.  Finally,  figures 
are  presented  to  show  the  percent  (%)  exceedance  vs.  percent  (%)  occurrence  at  each  bin  for 
significant  wave  height,  period,  and  direction,  respectively.  The  %  exceedance  does  not  always 
reach  100%  because  the  analysis  ignored  some  of  the  data  beyond  the  thresholds  applied  by  the 
user. 
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Figure  3.  Wave  Direction 

SUMMARY:  This  CHETN  describes  a  suite  of  tools  to  develop  wave  statistics  and  representative 
wave-climate  datasets  for  application  as  wave  forcing  at  offshore  boundaries  for  numerical  models. 
This  CHETN  is  the  first  of  a  two-part  series  discussing  the  binning  of  wave  conditions  for  use  in 
CMS-Wave  (or  other  spectral  wave  models  such  as  STWAVE).  The  second  CHETN  will  discuss 
importing  wave  conditions  into  CMS-Wave  and  using  CMS-Wave  output  as  input  wave  forcing  in 
GenCade.  Long-term  wave  information  from  hindcast  databases  or  measured  wave  information 
may  be  analyzed  using  the  tools  described  herein  to  develop  representative  wave  climate  statistics. 
These  results  are  applied  as  boundary  forcing  to  numerical  models  for  additional  wave 
transformation  and  as  input  to  predictive  tools  for  calculating  coastal  processes.  These  tools  may 
be  used  as  a  means  of  converting  measured  waves  and  model  solutions  to  manageable  model  input 
files.  These  tools  will  be  further  developed  and  integrated  with  existing  models  and  GUIs. 
Guidance  is  anticipated  to  change  following  publication;  the  user  is  directed  to  the  Coastal  Inlets 
Research  Program  (CIRP)  website  ( cirp. usace. army. m il)  and  the  CIRP  wiki  (cirp.nsace. 
armv.mil/wiki). 

POINTS  OF  CONTACT:  This  CHETN  was  prepared  as  part  of  the  Coastal  Inlets  Research 
Program  (CIRP)  and  was  written  by  Kenneth  J.  Connell  ( Kenneth  ConnellCwsolder.com)  of 
Golder  Associates,  Inc.,  and  Rusty  Permenter  ( Rusty. L. PermenteKmisace. army. mil  voice:  601- 
634-2089)  of  the  US  Army  Engineer  Research  and  Development  Center  (ERDC),  Coastal  and 
Hydraulics  Laboratory  (CHL).  Ashley  Frey  and  Tanya  M.  Beck  provided  peer-review  of  this 
publication.  For  information  about  CIRP,  please  contact  the  CIRP  Program  Manager,  Dr.  Julie 
Dean  Rosati  (tel:  202-761-1850)  or  by  email  (Julie. D.RosatiCcvusace. armv.mil).  This  CHETN 
should  be  referenced  as  follows: 
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NOTE:  The  contents  of  this  technical  note  are  not  to  be  used  for  advertising,  publication, 
or  promotional  purposes.  Citation  of  trade  names  does  not  constitute  an  official 
endorsement  or  approval  of  the  use  of  such  products. 
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